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Abstract 

We study the dynamics of nonlinear localized excitations ("soli- 
tons") in two-dimensional (2D) Bose-Einstein condensates (BECs) 
with repulsive interactions, loaded into an optical lattice (OL), which 
is combined with an external parabolic potential. First, we demon- 
strate analytically that a broad ("loosely bound", LB) soliton state, 
based on a 2D Bloch function near the edge of the Brillouin zone (BZ), 
has a negative effective mass (while the mass of a localized state is pos- 
itive near the BZ center). The negative-mass soliton cannot be held 
by the usual trap, but it is safely confined by an inverted parabolic 
potential ( anti-trap). Direct simulations demonstrate that the LB 
solitons (including the ones with intrinsic vorticity) are stable and 
can freely move on top of the OL. The frequency of elliptic motion of 
the LB-soliton's center in the anti-trapping potential is very close to 
the analytical prediction which treats the solition as a quasi-particle. 
In addition, the LB soliton of the vortex type features real rotation 
around its center. We also find an abrupt transition, which occurs 
with the increase of the number of atoms, from the negative- mass LB 
states to tightly bound (TB) solitons. An estimate demonstrates that, 
for the zero- vorticity states, the transition occurs when the number of 
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atoms attains a critical number N CT ~ 10 3 , while for the vortex the 
transition takes place at N CT ~ 5 • 10 3 atoms. The positive-mass LB 
states constructed near the BZ center (including vortices) can move 
freely too. The effects predicted for BECs also apply to optical spatial 
solitons in bulk photonic crystals. 

1 Introduction 

Periodic potentials, usually called optical lattices (OL), which are induced by 
interference of laser beams illuminating a Bose-Einstein condensate (BEC), 
are widely used in experimental studies of BECs. Not only one- dimensional 
(ID) OLs, but also their multidimensional counterparts are available in the 
experiment, see Refs. and references therein. 

An important type of nonlinear excitations that BECs may support are 
bright solitons (in this paper, following the trend commonly adopted in the 
current literature, we apply the term "soliton" to robust localized objects, 
without implying integrability of the underlying mathematical model). In an 
effectively ID ("cigar-shaped") condensate, solitons have been observed in 
the experiment [2]. Creation of solitons in a BEC confined in the OL has not 
been reported yet, although splitting of wave packets in a 87 Rb condensate, 
loaded in a cylindrical trap with an superimposed ID OL potential, which was 
reported recently j3] , might be a manifestation of transient soliton dynamics. 

Parallel to the experiments, many theoretical predictions for BEC soli- 
tons in OLs have been reported. In particular, attention was attracted to 
a possibility to create gap solitons, that are expected as a result of the in- 
terplay of the OL potential and repulsive interaction between atoms in the 
BEC jU Ej. While only ID solitons have been thus far observed in BECs, 
it was predicted that OLs may support multidimensional solitons of the gap 
type as well [4 j. A prediction which is especially relevant to feasible experi- 
ments is that OLs can also stabilize multidimensional solitons in BECs with 
attraction between atoms (because of the possibility of collapse, such solitons 
are unstable without a support potential) [6 . Moreover, 2D bright solitons 
carrying intrinsic vorticity can be stabilized too by the OL in the 2D case 
jHllZj. Various other properties of OL-supported solitons were studied too 


Similar localized objects are expected in a different physical setting, in 
the form of 2D optical spatial solitons in nonlinear photonic crystals (NPCs) 
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[HI El E] In optical media, an effective lattice potential is induced by periodic 
modulation of the local refractive index in the transverse direction(s) fOJ El 
Ej. An alternative is to use a virtual lattice, induced in a photorefractive 
medium by a transverse array of strong coherent laser beams; the latter 
arrangement has made it possible to observe 2D spatial solitons in a real 
experiment [TT], as well as their vorticity-carrying counterparts ^2] • 

Most theoretical studies of solitons in the above-mentioned models did not 
address their motion. In a recent work we considered in detail the mo- 
tion and collisions of ID solitons in the OL-equipped model in both cases of 
the attractive and repulsive cubic nonlinearity. In that work, the OL poten- 
tial was combined with an external parabolic (harmonic) trapping potential, 
which corresponds to the real situation in any BEC experiment. Direct sim- 
ulations of the corresponding effective ID Gross-Pitaevskii (GP) equation 
had demonstrated that, in either case (self-focusing or self-defocusing non- 
linearity), the soliton travels freely through the lattice if its norm (number 
of atoms) is below a certain critical value; above this value, it gets trapped 
by the OL. Depending on parameters, collisions between moving solitons are 
either elastic, or lead to their fusion. The most intriguing result concerns 
the motion of the soliton in the parabolic trap: in the case of attraction, the 
soliton behaves, as it may be naturally expected, as a quasi-particle in the 
harmonic-oscillator potential. However, in the case of repulsion, the soliton 
of the gap type is expelled by the parabolic trap, while it is stably confined 
by an anti-trap (repulsive parabolic potential). This result, which calls for 
direct experimental verification, was supported by analytical considerations 
which demonstrate that the effective mass of the gap soliton, regarded as a 
quasi-particle, is negative (in the case of the usual soliton, corresponding to 
the attractive nonlinearity, the effective mass is positive). 

An issue of straightforward interest is to extend these results to the 2D 
case, which is quite relevant, as experiments with 2D solitons in BECs seem 
possible right now. Besides that, the problem is also of direct relevance for 
steering spatial optical soliton in NPC media. In this work, we focus on 
the 2D BEC model with repulsive interactions, as this case is most relevant 
to the experiment, and produces most interesting physical results. First, we 
demonstrate that broad localized states, which we call loosely bound (LB) 
ones, may be represented as a product of a linear Bloch function, defined by 
the OL, which is taken close to the edge to the corresponding Brillouin zone 
(BZ), and a slowly varying envelope function, which obeys an asymptotic 
nonlinear equation that does not contain OL. Such a representation makes 
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it possible to demonstrate analytically that the effective mass of the LB 
state is negative. This way, we find not only fundamental LB states, but 
also ones carrying vorticity, whose mass is negative too; the existence of 2D 
localized vortices in the repulsive model has not been reported before. In 
direct simulations, the LB vortex features actual rotation around its center. 
Simulations confirm the stability of the localized LB states, and demonstrate 
that they can move freely through the OL. In particular, both fundamental 
and vortex LB states can be readily set in elliptic circular motion in the anti- 
trapping potential, the frequency of the circulations being very accurately 
predicted by the asymptotic equation. 

Direct simulations (as well as the asymptotic equation) demonstrate an 
abrupt transition (without hysteresis) from the LB states to much stronger 
localized tight-bound (TB) ones with the increase of the norm (number of 
atoms). The TB states are nothing else but the 2D gap solitons, that were 
first found in Ref. [3]. The LB-TB transition also takes place for the vortices, 
thus revealing the existence of new TB vortex solitons of the gap type. It is 
additionally found that the TB vortices undergo an intrinsic transition to a 
more symmetric shape with subsequent increase of the norm. Unlike their 
LB counterparts, the TB solitons cannot move, being strongly pinned by the 
OL. 

We also construct LB states, starting from the linear Bloch function close 
to the BZ center. In that case, the effective mass is positive, no transition 
to TB solitons occurs (as they do not exist), and the LB states, both fun- 
damental and vortex ones, are mobile, which is accurately described by the 
corresponding asymptotic equation. 

The rest of the paper is organized as follows. The analysis based on the 
derivation of the asymptotic equations for the LB states, with both negative 
and positive mass, is presented in section 2. Numerical results for static 
solitons (however, including rotation of the vortex LB states), such as the 
LB-TB transitions, are summarized in section 3. In section 4, we present 
numerical and analytical results concerning motion of the LB states, for the 
cases of the negative and positive mass. The paper is concluded by section 
5. 
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2 Formulation of the model and the analyti- 
cal approximation 

2.1 The Gross-Pitaevskii equations 

The mean-field description of the BEC dynamics is based on the GP equation 
for the mean-field wave function ip in three dimensions |14j . 



dip \ h 2 



_^V 2 + U (r) + G|Vf 
2m 



V>, (1) 



where m is the atomic mass, G = 4irh 2 a/m, a is the s -wave scattering length, 
U is the potential applied to the condensate, and the number of atoms is 
/ |^>| 2 dr. As is well known [3J E3 5 i n the case of a pancake-shaped con- 
figuration (when the frequency characterizing the confining potential in the 
transverse direction is much larger than the confining frequencies correspond- 
ing to the perpendicular directions), Eq. (JTJ) can be reduced to a 2D equation. 
After appropriate rescalings, such as t' = Ht / (m\ 2 ) , x' = x/X,y' — y/X and 
4> = G 1 / 2 ^, where A is the period of the spatially periodic potential, it takes 
the form 

i(fh = "^V 2 + ||0| 2 - e [cos(27nr) + cos(27n/)] + l -B{x 2 + y 2 )} 0, (2) 

where <fi is the renormalized wave function, the sign ' denoting the rescaled 
coordinates x',y',t' is omitted, the kinetic-energy operator V 2 acts on the 
coordinates x and y, the repulsive interaction between atoms is assumed, the 
OL period is normalized to be 1 in both directions (we consider the simplest 
case of the square lattice, although more sophisticated patterns, such as 
triangular, hexagonal, and quasi-periodic ones, may be interesting too), while 
e and B are amplitudes of the, respectively, OL and parabolic potential (in 
simulations, we will fix e = 10 and B = —0.005, or e = 5 and B = +0.005, 
negative B implying the anti-trapping potential). We will show the numerical 
results with the dimensionless equation (2). The unit length corresponds to 
the spatial period A of the optical lattice, and the unit time corresponds to 
m\ 2 /h ~ 1.3 x 10~ 3 s for A ~ 10~ 6 m and m( 87 Rb) ~ 1.37 x 10~ 25 kg. 

If B = and the nonlinear term is neglected, Eq. (j2J) becomes the 
ordinary 2D linear Schrodinger equation with the periodic potential, which 
has Bloch-state solutions. The Bloch function F(x, y) obeys the equation 

fiF(x, y) = ^V 2 F + e [cos(2vra;) + cos(27n/)] F, (3) 
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where /x is the chemical potential (or simply energy, in the case of a single 
atom), and the solution is quasi-periodic, i.e., F(x, y) = exp {ik x x + ik y y) f(x, 
where k x>y are two wavenumbers, and the function f(x, y) is periodic, such 
that f(x + 1, y) = f(x, y), f(x, y + 1) = f(x, y). The Bloch-wave solution de- 
termines the energy-band structure, \i = [J,(k x , k y ), the point (k x , k y ) = (ir, 7r) 
being the edge of the BZ (Brillouin zone). 

Following the lines of the analytical approach to the ID case developed in 
Ref. [T^j, an approximate solution to the full time-dependent GP equation 
((21) for a small- amplitude broad localized state ( "soliton" ) may be sought for, 
in the lowest approximation, as 

<i>{x,y,t) = e**F(x,y)$(x,y,t), (4) 

where F(x,y) is the Bloch function defined above, and <&(x,y,t) is an en- 
velope function which varies slowly (as a function of x and y) in compar- 
ison with F(x,y). This ansatz actually assumes that the nonlinearity and 
parabolic potential in Eq. (J2J) are small perturbations in comparison with 
the OL terms. The most typical gap soliton corresponds to a solution which 
approaches the linear Bloch function with k x = k y = it when its norm (the 
number of atoms), iV = f_™ \4>\ 2 dxdy, tends to zero (cf. a similar situation 
in the ID case [T5]). 

In the general case, the slowly varying amplitude $ obeys an asymptotic 
NLS equation, that can be derived by substituting the ansatz (JU) into Eq. 
(J2J) and applying a known averaging procedure 0j: 

Here, the effective mass M e g for linear excitations is determined from the 
curvature of the energy-band structure of the linear Bloch states, = 
(1/2) (d 2 /dk 2 x + d 2 /dk 2 y ) n, and 

^ Qti\F{x )y )\ A dxdy 

titi\F{x,y)\ 2 dxdy- {> 

The energy functional corresponding to the simplified GP equation (jHJ) is 

E = ^J J dxdy [M^ 1 1 V$| 2 + g |$| 4 + B (x 2 + y 2 ) |$| 2 ] . (7) 
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Stationary solutions to the asymptotic equation © can be looked for in 
an obvious form, 

$ = f(r) exp [i (mO — fit)] , (8) 

where r and 9 are polar coordinates in the plane (x,y), m = 0, 1,2, ... is a 
possible vorticity of the solution, and the amplitude f(r) obeys an equation 

f + r -if _ m 2 r - 2 = _ 2Meff ^ _ (jB/2) r 2 _ p/2 ] ; (g) 

[in the general case - in particular, in the case of the full GP equation ( |2J) 
which is not isotropic - the vorticity of a stationary solution may be defined 
as A9/ (27r), where A9 is the total change of the solution's phase along a 
closed path around the center]. The linear limit of Eq. Q is the linear 
Schrodinger equation for the 2D harmonic oscillator. For fixed integer m, 
the ground-state solution to the latter equation is (in the present notation) 

f(r) = r m exp(-aV/2), 



where a 2 = y/BM c &, with fi = (m + i)u, and uj = J B/M e s. 

2.2 Solutions near the edge and center of the Brillouin 
zone 

Near the BZ edge, i.e., close to the point (k x , k y ) = (tt,7t), the linear Bloch 
function can be approximated by a combination of four harmonics: 

F(x, y) = c\ exp(ik x x + ik y y) + C2 exp {i{k x — 2n)x + ik y y) 

+ C3 exp [ik x x + i(k y — 2ir)y) + C4 exp {i{k x — 2n)x + i{k y — 2ir)y) . (10) 

The substitution of this into Eq. (jSJ) and a usual truncation procedure (ne- 
glecting higher-order spatial harmonics) lead to an expression fi(k x , k y ) — 
/4 edge) + [(k x - vr) 2 + (k y - vr) 2 ] /(2M^ dse) ), where 

^ edge) = (vr 2 ±k|)/2, (11) 
In — \e\ 

and g = (3/4) 2 [recall g is defined by Eq. (JBJ)]; the signs ± in Eq. (|TT|) 
pertain to two Bloch-wave solutions which differ by a phase shift relative to 
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the periodic OL potential in Eq. J3J). As is seen, the effective mass in Eq. 
(|12|) is negative for the linear excitations near the BZ edge (unless \e\ is very 
large, in which case the above approximation does not apply). 

Equation © with M e s < 0, g > 0, and B < (anti-trap) is tantamount to 
the usual 2D GP equation with positive mass, negative g (which corresponds 
to self-attraction), and B > (usual trap, rather than anti-trap). As is well 
known fUJ E| j the GP equation in the latter form has localized solutions, 
unless the norm iV is too large (otherwise, collapse will take place). These 
solutions, and their comparison with numerically found solutions to the full 
underlying equation (J2J), will be displayed below. 

An approximation similar to that based on Eq. (|1U|) can also be developed 
close to the BZ center. At the BZ center, the linear Bloch function may be 
approximated by 

F(x,y) = [l + c 2 cos(2ttx)] [1 + c 2 cos(2rty)] , (13) 

and the eventual expressions for the effective mass and /io are 

^center) = 2tT 4 + S 2 + 71 V4tT 4 + 2e 2 

eff ~ lOvr 4 + e 2 - 3n 2 V^ 4 + 2e 2 ' 1 ' 



^(center) = ^ _ ^4 + ^2/2. ( i5 ) 

Other coefficients are found to be c 2 = —^ofs and g = (1 + 12c! + 6c 4 ) /(l + 
2c|). In this case, the effective mass is always positive [in particular, in 
the case with e = 5 and B = 0.005, for which this approximation will be 
compared to numerical results below, M e g 6nter ^ pa 1.12]. Obviously, Eq. Q 
with M e ff > 0, g > 0, and B > has a localized solution (actually, of the 
Thomas-Fermi type, see below) for any vorticity m and norm N. 



3 Numerical results 
3.1 Loosely bound states 

Our first objective is to find LB ("loosely bound") localized states near the 
BZ edge, which can be constructed using the asymptotic GP equation (0) 
and the relations (@J) and (fiT)|) . and compare them with the corresponding 
numerical solutions to the full equation (j2J). The 2D localized solutions to 
the asymptotic GP equation (j3J) were found numerically f0JE3, and in an 
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approximate form they can be also solved by dint of a variational method 
[T7j . In particular, the radial equation (0) can be solved numerically by 
the shooting method, which simultaneously makes it possible to find the 
nonlinear eigenvalue \x. Figure 1 displays numerical solutions to Eq. (0 ) for 
m — 0, 1 and 2, obtained by means of the shooting method. The norm of the 
corresponding solutions for the underlying wave function <p is, according to 
Eq. N m J J f 2 (r) cos 2 (irx) cos 2 (7ry)dxdy ~ N /4, where N is the norm 
of $ [see Eq. ©]. In particular, N = 1.52 for m = 0, iV = 2.98 for m = 1 , 
and N = 4.08 for m = 2. 




Figure 1: Typical examples of numerically found solutions to Eq. with 
M eff = -10/(2tt 2 - 10), g = (3/4) 2 and B = -0.005 which represent 2D 
solitons with vorticity m — 0, 1 and 2. 

To characterize families of the solutions, Figs. 2(a) and 2(b) display, 
respectively, the absolute value of the energy-per-atom, — e = —E/N for the 
solution [the energy E is given by the functional (7)], and the maximum 
value of /(r) vs. the norm N. In particular, — e decreases with N, because 
a larger atom density gives rise a larger positive contribution to E through 
the term g\&\ 4 in the energy density. The maximum value of f(r) increases 
with N, since the effective interaction is attractive owing to the negativeness 
of M eff . 

Figures 3(a), (b) and (c) present comparison of the above semi-analytical 
solutions (for m — 0, 1 and 2), based on Eqs. (jU), (fTUj) and (J3J, with di- 
rect numerical solutions of the full 2D equation ( |2J). The comparison is 
performed by displaying the semi-analytical and direct numerical solutions 
for \<f)(x, L/2)\, taken along the central cross-section, at y = L/2 (the norms 
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Figure 2: The absolute value of the energy-per-atom, — e, (a) and amplitude 
(b) of the 2D solitons with the vorticity 0, 1 and 2, as found from Eq. (jUJ), 
vs. the norm. 

of the numerical and semi-analytical solutions are virtually identical in each 
case). Fairly good agreement is seen for m — 0, 1 and 2, although \<fi\ in the 
direct numerical solutions is slightly larger. 

From the viewpoint of the reduced GP equation (J3j), which does not 
explicitly contain the OL potential, the angular velocity Q = fi/m, which 
appears in the expression (JHJ), is only a phase velocity. However, the relation 
(J3J) suggests that the corresponding solution to the full GP equation (J2J) may 
feature actual rotation of the pattern with the angular velocity Q. Direct 
simulations confirm this expectation, as shown in Fig. 0J The figure displays 
three consecutive snapshots of the region in which the vortex solution with 
m = 2 shown in Fig. 3(c) takes sufficiently large values, namely, satisfying 
conditions [cos(7ra;) cos(iiy)} Recj)(x, y) > and \<p(x, y)\ > 0.06. The time 
interval between the snapshots is 1. This figure clearly shows that the vortex 
state indeed rotates with the frequency Q m fi/m in the OL potential. This 
figure also demonstrates that the vortex state with m = 2 is stable. The 
vortex solutions, with m = 1 and 2, to the asymptotic equation (jSJ) are known 
to exist PH] and be dynamically stable [HI if N is smaller than a critical value. 
We have also checked the stability of the vortex states with m = 1 and 2, 
shown in Fig. 3, by direct numerical simulations of Eq. (2), starting from 
initial conditions in the form of perturbed vortex states. When the norm 
iV is increased, more complicated behaviors are expected. Some simulation 
results such as a transition to a tightly bound state and a symmetry-breaking 
process are shown in the next section. 

Solutions predicted by the asymptotic equation (jSJ) near the BZ center, 
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Figure 3: The central cross section of solutions \<f>\ for the loosely bound 2D 
states near the edge of the Brillouin zone with (a) m = 0, (b) m = 1, and 
(c) m = 2, which corresponds to Eq. (fTUj) . Solid and dashed curves depict, 
respectively, the solutions obtained from the asymptotic equation (J3J) and 
full GP equation ©. 

i.e., corresponding to Eqs. (JTHJ) and (JHJ), were also found, and compared to 
direct numerical solutions of the underlying GP equation (j2J). It should be 
stressed that, as the effective mass (|14|) is positive in this case, the correspond- 
ing localized solutions to the asymptotic equation with g > (repulsive 
interaction) are not soliton-like ones, but are instead close to the Thomas- 
Fermi (TF) states, that can be obtained neglecting the kinetic-energy term 
in the GP equation [Hj. Figure El shows typical examples of thus found so- 
lutions for m — 0, m — 1, and m — 2. As is seen, the accuracy provided by 
Eqs. (0) and (j!3|) is very good in this case, for all the values of m. 

3.2 Relation to tightly bound solitons 

The full GP equation (J2J cannot be reduced to the asymptotic NLS equation 
(J3J) if the assumption of the slow variation of the envelope function is not 
valid; in that case, strongly confined ("tightly bound", TB) states may be 
possible, that are not described by Eq. (J3J). To investigate such solutions, 
and compare them to their loosely bound counterparts, we performed direct 
numerical integration of Eq. (J2J), looking for solutions that could be classified 
as having the vorticity m = or m = 1 [in terms of Eq. (jHj)]. Figures EJa) 
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(a) (b) (c) 




Figure 4: An example of the uniform rotation of the vortex state with N = 
4.08 and m = 2. 

and|H^b) display the amplitude of the thus found solutions [the maximum of 
\4>(x, y)\] as a function of the norm, N = J J \(p\ 2 dxdy. The lower branches 
in these figures precisely correspond to the LB states described above, the 
dashed curves showing the dependence of the amplitude vs. N as found 
from the asymptotic equation (jSJ). The solutions which belong to the upper 
branches in Fig. 6 can be identified as TB (strongly localized) 2D gap-soliton 
states in the GP equation with repulsive interaction between atoms, that 
persist in the limit of B = and were first found in Refs. jlj. There is weak 
hysteresis for m = between the LB and TB states for the parameter values 
of e = 10 and B = —0.005. The discontinuous transition between the LB 
and TB states occurs for sufficiently large e. (The discontinuous transition 
was not observed for e=5 and B = —0.005.) Figures 6(a) and (b) also make 
it obvious that the asymptotic equation (j3J) predicts a jump upward from the 
lower branch, which corresponds to the collapse in the asymptotic 2D NLS 
equation ©. However, the collapse does not occur in the full GP equation 
((21) with the repulsive nonlinearity (note that the negative effective mass is 
irrelevant when the localization length is on the same order of magnitude 
as the OL period). The critical value of the norm for the collapse in the 
asymptotic equation (jHJ) is somewhat larger than the value at which the LB- 
TB transition occurs in the numerical solution of Eq. (j2*j). 

It is necessary to estimate the actual number of atoms in the BEC at 
which the LB-TB transition, shown in Fig. 6, is expected. Assuming a "pizza- 
shaped" condensate, with the diameter ~ 1000 OL periods and transverse 
width ten times smaller, taking the scattering length 5 nm (which is close to 
its value for 87 Rb), and following once again the derivation of the 2D equation 
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p (a) p (b) p (c) 




Figure 5: Dashed curves correspond to Thomas-Fermi-like solutions of the 
asymptotic equation (jSJ) with B = +0.005 (the usual trap, rather than anti- 
trap) and positive effective mass, M e ^ = 1.12, as found by the shooting 
method from Eq. (9). These solutions describe loosely bound states close 
to the center of the Brillouin zone, as per Eq. ( [TSJ). Continuous curves are 
numerical counterparts of these solutions, obtained from the full GP equation 
(12) with e = 5. (a) m = 0, N = 19.7; (b) m = 1, JV = 29.8; (c) m = 2 , 
N = 33.0. 

(J2J) from the full 3D GP equation, we conclude that the real number of atoms 
is obtained by multiplying the numbers on the horizontal axes in Fig. 6 by a 
factor ~ 10 3 . 

To further illustrate the drastic difference between the LB and TB states 
with zero vorticity (m = 0), we display their typical examples, found as 
direct numerical solutions of Eq. (J2J), in Figs. 7(a) and (b). These figures 
display regions in which the absolute value of the respective solution exceeds 
the quarter- amplitude, i.e., \<f>(x, y)\ > |0| ma x/4. Additionally, Fig. 7(c) 
displays the absolute value for both solutions, |0( function of x 

along the diagonal, x = y. 

In a similar manner, the comparison of the LB and TB states with 
m — 1 is presented in Fig. 8, where the panels (a) and (b) again display 
regions where, respectively, the loose and tight solutions take the values with 
\4>(x, y)\ > |0| max /4, and the panel (c) shows the diagonal cross sections of 
both solutions. In particular, the TB vortex in Fig. 8(b) features a four-fold 
symmetry (it is invariant with respect to the rotation by vr/2). In this con- 
nection, it is relevant to mention that, while the same symmetry is observed 
in all the TB vortices for a sufficiently large norm, N > 10, the four-fold 
symmetry is broken for N < 9. 

To illustrate the latter property, Fig. 9 displays three consecutive snap- 
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Figure 6: The amplitude of the 2D solution with m — (a) and m = 1 
(b) vs. the norm (number of atoms) N, in the case of e — 10 and B = 
—0.005 [the dashed lines show the same dependencies, as predicted by the 
asymptotic equation © for the loosely bound states]. The lower and upper 
branches of the dependence correspond to the loosely and tightly bound 
states, respectively (see further details in the text). 

shots of regions with \<f>\ > 0.5, taken from the simulation starting with an 
LB vortex state that has m = 1 and N = 7. As is seen, the patterns do 
not feature the four-fold symmetry, unlike Fig. 8(b), but a nearly two-fold 
symmetric TB state is generated. 

The transition between the less symmetric and more symmetric TB vor- 
tices accounts for the slope breaking in the amplitude-vs.-norm curve for 
N = 9, which is evident in Fig. 6(b). We have also found the TB vortices 
with m = 2 and checked that they are stable as well. Figure 10 displays 
the TB vortices with m = 2 for N = 3.92 and N = 48.1. Figures 10(a) 
and (b) display the regions where the TB solutions take the values with 
\4>(x, y)\ > |0|max/4, and the panel (c) shows the cross sections along the 
line y — L/2 = — 0.45 (a; — L/2), on which |</>| takes the maximum value for 
N = 48.1. To the best of our knowledge, the existence of TB vortex solitons 
in the 2D model with self-repulsion has not been reported before, there- 
fore theses result (displayed here for m = 1 and 2) are novel by themselves 
too (similar stable vortices were very recently found by means of a different 
method HH). 

The localized solutions belonging to the LB and TB branches differ dras- 
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Figure 7: Typical examples of the loosely (a) and tightly (b) bound 2D states 
with m = (zero vorticity) for e = 10 and B = —0.005. The norms of the 
two states are, respectively, Moose = 0.393 and iVtight = 2.66. Panel (c) 
displays the diagonal cross sections of the two solutions, \<fi(x, x)\ (the solid 
and dashed lines depict the loosely and tightly bound states, respectively). 

tically not only in the dependence of |0| m ax vs. N (i.e., in their static prop- 
erties), but in dynamical properties too. As it will be shown in the next 
section, the TB solitons are always pinned by the underlying OL potential, 
while the LB states may move freely. In fact, a sharp transition between 
mobile and immobile ID solitons with M c r < (i.e., close to the BZ edge) 
and g > (repulsive interactions) was reported in Ref. ^Hj- As well as in 
the 2D case, the immobilization takes place with the increase of the soliton's 
norm, and no hysteresis between the mobile and immobile solitons is found. 

The abrupt transition between the LB and TB localized states, as de- 
scribed above, is a characteristic feature of the localized states with the 
negative effective mass, found close to the BZ edge, where the LB state is 
well approximated by Eqs. (jlj), ©, an d (fTU|) . On the other hand, no simi- 
lar transition is found for the localized states close to the BZ center, where 
the approximation (fTU|) is replaced by (|T3*|) . The absence of the distinction 
between LB and TB states in the latter case is explained by the fact that, 
as is commonly known, TB solitons simply do not exist in the case of the 
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Figure 8: The same as in Fig. for the loosely (a) and tightly (b) bound 
vortices with m — 1. The norms of the two states are, respectively, N\ oose = 
0.413 and N tight = 13. 

self-repulsion (g > 0) near the BZ center. Localized LB states of course exist, 
as suggested by the TF approximation. 



4 Motion of loosely bound localized states 

In the 2D model with both repulsive [3] and attractive [H] interaction, the TB 
solitons are strongly pinned by the OL potential, and cannot move (stable 
TB solitons also exist in the 2D attractive model with a quasi-lD periodic 

(a) (b) (c) 
10 j— , 1 i i (0 j— j j 10 i 1 , , 
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Figure 9: Three snapshot patterns taken at the moments (a) t = 10, (b) 
t = 40, and (c) t = 80 in the simulation of a loosely-bound vortex with 
m = l, evolving into a tightly-bound state. The norm is N = 7. 
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Figure 10: The same as in Figs. [7|and|B]for the loosely (a) and tightly (b) 
bound vortices with m = 2. The norms of the two states are, respectively, 
Moose = 3.92 and Alight = 48.1. (c) A profile of \<p(x, y)\ for the tightly-bound 
vortex along the line y — L/2 = — 0A5(x — L/2). 

potential, and in the 3D model with a quasi-2D potential, in which case 
they may freely move in the unconfined direction [20J). On the contrary 
to that, the LB localized states described above may move freely on top of 
the underlying 2D periodic potential, keeping their coherent structure. This 
feature resembles a similar one that we have recently demonstrated in the 
ID version of the present model ^3]; i n particular, the ID solitons with the 
negative effective mass are expelled by the ordinary parabolic trap, and are 
held by the anti-trap. However, the character of the motion in the 2D case 
may be completely different, see below. 

A 2D LB state can be set in motion either by an initial shift from the 
central position [relative to the (anti-)trap], or by lending it an initial mo- 
mentum. We have performed numerical simulations, using both ways to push 
the localized state. For instance, Fig. 10 displays results obtained with the 
initial condition <p(x, y, 0) = 4>o(x — 5, y) exp (ikq y y), where 4>o(x, y) is the sta- 
tionary LB state with m = and A" = 0.479, whose peak position is located 
at the central point, (L/2, L/2), and q y = 0.2. In other words, this initial 
condition implies that the localized state is pushed by shifting its center in 
the x direction, and giving it an initial momentum along y. 

We define the instantaneous position of the center of the localized state as 
X = // x|0(x, y; t)\ 2 dxdy/ J J \<p\ 2 dxdy. Figure 11(a) demonstrates that, as a 
result, the center of the LB state moves along an elliptic trajectory [in the case 
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displayed in Fig. 11(a), the simulations were run up to the time t max = 100, 
while the numerically measured period of the elliptic motion was T num pa 88, 
i.e., the localized state has completed one full circulation (in the clockwise 
direction, since the initial velocity v y = q y /M e s is negative). To make the 
difference from the strongly pinned TB solitons clearer, Fig. 11(b) displays 
a result of a similar attempt to set in motion the localized state of the TB 
type. As a result, initially its center moves irregularly, and is then trapped 
by the OL potential. Many other runs of the simulations produced quite 
similar results; in particular, the eccentricity of the elliptic orbit depends 
on the initial push, and effectively ID periodic motion along straight lines 
through the center of the domain are possible too. In this connection, it is 
relevant to mention that, although the underlying OL makes the medium 
anisotropic, simulations do not reveal any dependence of the frequency on 
the direction of the ID oscillations (the frequency of the elliptic motion does 
not depend on the orientation of the ellipse either). 




Figure 11: Panel (a) displays a typical example of the stable elliptic motion 
of a localized loosely-bound state with zero vorticity (m = 0) and negative 
effective mass (M eff pa —1.03) in the anti-trap with B = —0.005. In this case, 
e = 10, and the norm of the localized state is N = 0.479. For comparison, 
panel (b) displays an attempt to set in motion a tightly-bound soliton with 
N = 4.01. 

The periodic motion of the LB state can be easily described within the 
framework of the asymptotic equation (j3J). In that case, one can demonstrate 
that motion for the center X of the localized state, which is treated as a rigid 
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quasi-particle, obeys the Newton's equation of motion, 



M eff 



d 2 X 



-BX 



(16) 



dt 2 



[M e ff is the same effective mass as in Eq. (jSj), while the right-hand side is the 
potential force generated by the parabolic (anti-)trap]. Equation (JTfiJl indeed 
describes the observed motion with a good accuracy. For instance, in the 
case displayed in Fig. El (a), Eq. (fi7)J) predicts the period of the motion, 



which should be compared to the above-mentioned numerically found period, 
T num « 88. 

As it was explained above, the localized states found close to the BZ center 
are actually always of the LB type, hence one may expect that they may 
move freely too. This is indeed observed in simulations. Typical examples 
are displayed in Fig. 12, for the localized states with m = and 2 (quite a 
similar motion mode was observed for the vortex with m = 1). The trajectory 
corresponds to anti-clockwise elliptic circulation in this the initial 

velocity v y = q y /M e s is positive. 

We stress that the vortex LB states readily exhibit circular motion with 
the same period as their counterparts with m — 0, and the analytical predic- 
tion for the circulation period, based on Eq. (JT7J), is very accurate in this case 
too. For example, in both cases (m = and m = 2) shown in Fig. 12, the 
predicted period is T ana i rs 94.0, while the one extracted from the numerical 
data is T num 94.5. 

5 Conclusion 

In this work, we have developed analysis of dynamics of nonlinear localized 
excitations in two-dimensional (2D) Bose-Einstein condensates (BECs) with 
repulsive interactions, which are subjected to the action of the optical lattice 
(OL) combined with the external parabolic potential. First, we have shown, 
in an analytical form, that broad ("loosely bound", LB) localized states, 
based on the 2D Bloch function near the edge of the Brillouin zone (BZ), 
feature a negative effective mass (while the mass of the LB localized states is 
positive near the BZ center). The negative-mass pulse cannot be confined by 




(17) 
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Figure 12: Panel (a) has the same meaning as in Fig. 11, but for the case 
of the positive effective mass, M eff pa 1.12, and in the ordinary trap with 
B = +0.005 (and e — 5). Panel (b) shows the corresponding time dependence 
of the central coordinate X of the moving localized state. In panel (c), the 
same dependence is shown for the circular motion of the vortex with m — 2. 

the usual trap, but it is held by an inverted parabolic potential (anti-trap ). 
Direct simulations have demonstrated that the LB states (including vortices) 
are stable and can freely move on top of the OL, and the frequency of their 
periodic motion along an elliptic orbit in the anti-trapping potential agrees 
very well with the analytical prediction, which treats the LB solitons as quasi- 
particles. Besides that, the LB vortices feature rotation around the center. 
With the increase of the number of atoms, an abrupt transition occurs, from 
the negative-mass LB states to their tightly bound (TB) counterparts; unlike 
the LB states, the TB solitons are rigidly pinned by the OL. The positive- 
mass LB states (including vortices), constructed near the BZ center, are 
stable and move freely too. 

Further dynamical effects may be considered in the 2D situation, such 
as collisions between moving solitons. In fact, various types of collisions 
are possible, such as between zero-vorticity solitons and vortices, as well as 
between LB and TB states with the negative mass (the latter is possible 
in the presence of the anti-trap). On the other hand, interactions between 
solitons with positive and negative mass are precluded, as only one type of 
the solitons may be confined by a given trapping or anti-trapping potential. 
The consideration of collisions is, however, beyond the scope of this paper. 
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We expect that similar LB and TB states may be found in BECs subject 
to the action of the 3D (three-dimensional) OL potential. However, simula- 
tions of the corresponding 3D GP equation is a rather difficult problem. 

Estimates for physical parameters at which the predicted effects may 
be observed are similar to those in the ID counterpart of the model [32J 
. In particular, a circular motion of the LB soliton in the anti-trap may 
be expected at a frequency ~ 100 Hz in a domain of the diameter ~ 0.5 
mm (~ 1000 spatial periods of the OL). The transition from the LB to TB 
solitons is expected when the number of atoms in the condensate attains the 
values between ~ 2 • 10 3 for the m = states and ~ 6 • 10 3 for the vortex 
with m = 1. 

In the experiment, the LB and TB localized states are expected to form 
spontaneously from a properly chosen number of atoms loaded in the appro- 
priate potential, as there is virtually no overlap between these two types of 
the states. As concerns vortex states, they may be created using the well- 
known technique of optical stirring. The effects predicted in this work for 
the BEC can also take place with spatial optical solitons in bulk photonic 
crystals. 
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